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Ab initio calculations of the quasi-elastic electromagnetic and neutral-weak response functions of 
4 He and 12 C are carried out for the first time. They are based on a realistic approach to nuclear 
dynamics, in which the strong interactions are described by two- and three-nucleon potentials and the 
electroweak interactions with external fields include one- and two-body terms. The Green’s function 
Monte Carlo method is used to calculate directly the Laplace transforms of the response functions, 
and maximum-entropy techniques are employed to invert the resulting imaginary-time correlation 
functions with associated statistical errors. The theoretical results, confirmed by experiment in 
the electromagnetic case, show that two-body currents generate excess transverse strength from 
threshold to the quasi-elastic to the dip region and beyond. These findings challenge the conventional 
picture of quasi-elastic inclusive scattering as being largely dominated by single-nucleon knockout 
processes. 

PACS numbers: 21.60.De, 25.30.Pt 


In first-order perturbation theory, the interactions 
of an external electroweak probe with a nucleus are 
described by response functions. These response 
functions—two for the processes A(e, e') induced by 
electromagnetic interactions, and five for the processes 
A{yi,u[) and A(Pi,u(), or A[yi,l~) and A(pi,l + ), in¬ 
duced by neutral or charge-changing weak interactions— 
determine the inclusive differential cross sections |T]. 
They can be written schematically as 

R a 0(q,u)~^26(uj + e q - E f ){f\0 a (q)\0)*{f\Op(q)\0), 
f 

(!) 

where q and oj are the momentum and energy transfers 
injected by the external field into the nucleus, |0) and |/) 
represent respectively its initial ground state of energy 
Eq and final continuum state of energy Ef, O a (q) and 
Oft (q) denote appropriate components of the of the nu¬ 
clear electroweak current operator (their w-dependence 
is dealt with as described below), and an average over 
the ground-state spin projections is understood (pre¬ 
cise definitions for the nuclear electroweak response func¬ 
tions, and resulting inclusive cross sections, are given in 
Ref. P]). 

At large values of momentum and energy transfers 
(q > 1 GeV and u > 0.5 GeV), where the dynamics of 
interacting nucleons is inextricably interwoven with the 
internal dynamics of individual nucleons, the accurate 
calculation of the response functions poses formidable 
challenges, particularly in view of the fact that a con¬ 
sistent theoretical framework to describe such a regime 
is still lacking. Even at the lower q and u of interest 
in the present study (q < 0.5 GeV and ui in the quasi¬ 
elastic region), where the consequences of the nucleon’s 
substructure on nuclear dynamics can be subsumed into 


effective many-body potentials and currents, this calcula¬ 
tion remains extremely difficult: it requires knowledge of 
the whole excitation spectrum of the nucleus and inclu¬ 
sion in the electroweak currents of one- and many-body 
terms. 

In the case of inclusive weak scattering, a further dif¬ 
ficulty exists for comparing calculated and experimen¬ 
tal results. The experimental initial neutrino energy is 
not known; instead there is a broad spectrum of ener¬ 
gies. This means that the observed cross section for a 
given energy and angle of the final lepton follows from 
a folding with the energy distribution of the incoming 
neutrino flux and, consequently, may include contribu¬ 
tions from q-co regions where different mechanisms are 
at play: the threshold region, where the structure of the 
low-lying energy spectrum and collective effects are im¬ 
portant; the quasi-elastic region, which is dominated by 
scattering off individual nucleons and nucleon pairs (see 
below); and the A-resonance region, where one or more 
pions are produced in the final state [2j. 

Integral properties of the response functions can be 
studied by means of sum rules, which are obtained from 
ground-state expectation values of appropriate combi¬ 
nations of the current operators (and commutators of 
these combinations with the Hamiltonian in the case 
of energy-weighted sum rules), thus avoiding the need 
for computing the nuclear excitation spectrum. Ab ini¬ 
tio quantum Monte Carlo (QMC) calculations of (non- 
energy-weighted) electroweak sum rules in 12 C have been 
recently reported in Refs. m ®- These calculations 
have demonstrated that a large fraction (~ 30%) of the 
strength in the response arises from processes involving 
two-body currents, and that interference effects between 
the matrix elements of one- and two-body currents play a 
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major role [5]. These effects are typically only partially, 
or approximately, accounted for in existing perturbative 
or mean-field studies m- 

Yet, sum rules do not provide direct information on 
the distribution of strength, whether, for example, the 
calculated excess strength induced by two-body currents 
is mostly at lar ge ui, we ll beyond the quasi-elastic peak 
energy to qe = q 2 + ?n 2 — m (m is the nucleon mass), or 
is also found in the quasi-elastic region with ui < w qe . 
Moreover, in the electromagnetic case, comparison of 
theoretical and experimental sum rules is problematic, 
since longitudinal and transverse response functions ob¬ 
tained from Rosenbluth separation of measured inclusive 
(e, e') cross sections are only available in the space-like 
region (w < q) and therefore must be extrapolated into 
the unobserved time-like region (w > q) before “experi¬ 
mental” values for the sum rules can be determined, see 
Refs, d nm for a discussion of these issues. 

In this paper we report the first ab initio calculations of 
the electromagnetic and neutral-weak response functions 
of 4 He and 12 C (other studies for 4 He have been already 
performed within different frameworks, i.e. m my 
These calculations proceed in two steps: the first involves 
the use of QMC methods to compute the response in 
imaginary time, the so-called Euclidean response mm, 
while the second consists in the inversion, via maximum 
entropy techniques mm, of these “noisy” imaginary- 
time data to obtain R a p(q,u>). The dynamical frame¬ 
work is based on a realistic Hamiltonian, including the 
Argonne uis two-nucleon 117] (AV18) and Illinois-7 three- 
nucleon [TS| (IL7) potentials, and on realistic electroweak 
currents with one- and two-body terms. A concise de¬ 
scription of this framework is in Refs. 13 m> while a 
more extended one can be found in the reviews da HO]. 
These latter papers also illustrate the level of quantita¬ 
tive success it has achieved in accurately predicting many 
properties of s- and p-shell nuclei up to 12 C, including, 
among others, energy spectra of low-lying states, static 
properties like charge radii, magnetic dipole and elec¬ 
tric quadrupole moments, radiative and weak transition 
rates, and elastic and inelastic electromagnetic form fac¬ 
tors. 

The Euclidean response function is defined as the 
Laplace transform 

pOO 

E a p(q,r) = C a p{q) duj e~ TU1 R a p(q,Lo) , (2) 

Ju th 

where w t h is the inelastic threshold and the C a p are q- 
dependent normalization factors. In R a p(q,oj) the in¬ 
dependence enters via the energy-conserving (5-function 
and the dependence on the four-momentum transfer 
Q 2 = q 2 — to 2 of the electroweak form factors of the 
nucleon and iV-to-A transition in the currents. We 
remove the latter by evaluating these form factors at 
Qqe = q 2 — w qe - In the case of the electromagnetic lon¬ 
gitudinal (L or a/3 = 00) and transverse (T or a/3 = xx) 


response functions, the normalization factors are m 
Cl = Ct = 1/ [G^(Qq e )] 2 , where G P E is the proton 
electric form factor, while in the neutral-weak response 
functions they are the same as those adopted in the sum 
rule calculations reported in Ref. [4]. With these defini¬ 
tions the response functions in Eq. ©> can be thought of 
as being due to point-like, but strongly interacting, nu¬ 
cleons. Note that non-energy-weighted sum rules 13® 
correspond to E a p(q 1 T = 0), while energy-weighted ones 
are obtained by taking derivatives of E a p(q,r) with re¬ 
spect to r and evaluating them at r = 0. 

The Euclidean response can be expressed as a ground- 
state expectation value, 

E a p(q,T) = (0|Ot(q) e -(g-SoA O/3 (q)|0) 

C a0 (q) (0\e-( H ~ E o)r\0) ’ 1 J 

where H is the nuclear Hamiltonian (here, the AV18/IL7 
model), r is the imaginary-time, and Eq is a trial energy 
to control the normalization. In this paper we report 
responses computed with the variational wave function, 
|0> = |\fv); in Refs. [3 EO it was shown that sum rules 
computed with |\Ev) are very close to those computed 
with the exact Green’s function Monte Carlo (GFMC) 
wave functions. The calculation of the matrix element 
above is carried out with GFMC methods [T3] similar to 
those used in projecting out the exact ground state of H 
from a trial state [21]. It proceeds in two steps. First, an 
unconstrained imaginary-time propagation of the varia¬ 
tional Monte Carlo (VMC) state |'lv) is performed and 
saved. Next, the states O j a(q)|’Iv) are evolved in imag¬ 
inary time following the path previously saved. During 
this latter imaginary-time evolution, scalar products of 
exp [— ( H — E 0 ) Ti} 0 / 3(q)|\I>v’) with O a (q)|\fv) are eval¬ 
uated on a grid of r ; ; values, and from these scalar prod¬ 
ucts estimates for E a p(q,Ti) are obtained (a complete 
discussion of the methods is in Refs. mm)- 

In Fig. [I] the electromagnetic longitudinal (El, top 
panel) and transverse (Et, lower panel) Euclidean re¬ 
sponse functions of 12 C are compared to those extracted 
from the world data analysis by Jourdan |22j , represented 
by the shaded bands. In order to better show the large r 
behavior, all the figures in this paper show E a p(q ) r) = 
exp[r q 2 /(2m)\E a p(q,T); this scaled response would be 
a constant for an isolated proton. The “experimental” 
E L (q,r) and ET(q,x) follow from Laplace-transforming 
the longitudinal and transverse data. These are first di¬ 
vided by \G p e (Q 2 )\ to obtain corresponding response 
functions of point-like nucleons, and then integrated with 
the weight factor exp(—rw) up to w max , where measure¬ 
ments are available. The strength in the unobserved re¬ 
gion with ui > Wmax is estimated by assuming that the 
R L (q,u > w max ) and R T (q,u > w max ) of 12 C are pro¬ 
portional to those in the deuteron, which can be ac¬ 
curately calculated [T]. The procedure is identical to 
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FIG. 1. (Color online) Euclidean electromagnetic longitudinal 
(top panel) and transverse (lower panel) response function of 
12 C at q = 570 MeV. Experimental data are from Ref. }22j . 

that used in Ref. 0 for the sum rules. As discussed 
in Ref. [3j, the scaling assumption can be justified by ob¬ 
serving that the high uj (well beyond w qe ) region of the 
response is dominated by two-nucleon physics, in partic¬ 
ular by deuteron-like np pairs in the ground-state of the 
nucleus. It is important to stress that, as r increases, 
the Euclidean response functions become more and more 
sensitive to strength in the quasi-elastic and threshold 
regions of Rl,t{q, w). Indeed, in this limit (r > l/w qe ) 
contributions from unmeasured strength at w > w max are 
exponentially suppressed. 

In Fig. [T] we show results obtained by including only 
one-body (open circles) or both one- and two-body (solid 
circles) terms in the electromagnetic transition operators. 
In the longitudinal case, destructive interference between 
the matrix elements of the one- and two-body charge op¬ 
erators reduces, albeit slightly, the one-body response. 
In the transverse case, on the other hand, two-body cur¬ 
rent contributions substantially increase the one-body re¬ 
sponse. This enhancement is effective over the whole 
imaginary-time region we have considered, with the im¬ 
plication that excess transverse strength is generated by 
two-body currents not only at w ^ w qe , but also in the 
quasi-elastic and threshold regions of Rr(q,oj). It is re¬ 
assuring to see that the full predictions for both longitu¬ 


dinal and transverse Euclidean response functions are in 
excellent agreement with data. 

At larger values of r the statistical errors associated 
with the GFMC evolution are rather large, particularly 
in the longitudinal response for which the elastic contri¬ 
bution proportional to the square of the 12 C form fac¬ 
tor [5] needs to be removed in order to account for the 
inelastic strength only. However, it should be possible 
to reduce these errors in the future by investing substan¬ 
tial additional computational resources in this type of 
calculation. Those presented here were performed with 
~45 million core hours of Argonne National Laboratory’s 
IBM Blue Gene/Q (Mira) parallel supercomputer. The 
Automatic Dynamic Load Balancing (ADLB) library [25] 
was used to distribute the imaginary time propagation of 
Ofl(q)ITv') and the evaluation of the matrix element in 
Eq. |3| over more than 8000 MPI ranks. The code is at 
present approximately 75% efficient at this scale. 

In Fig. [2] we show the largest of the five Euclidean 
neutral-weak response functions: the transverse (top 
panel) and interference (lower panel) E a p(q,r), having 
respectively a/3 = xx and a/3 = xy in the notation of 
Ref. [Tj. The E xy (q 1 r) response is due to interference 
between the vector (VNC) and axial (ANC) parts of the 
neutral current (NC), and in the inclusive cross section 
the corresponding R xy (q, u>) enters with opposite sign de¬ 
pending on whether the process A(yi,v[) or A(Vi,v[) is 
considered [QQ. On the other hand, in the transverse 
case the interference of VNC and ANC terms vanishes, 
and E xx (q,r) is simply given by the sum of the terms 
with both O a and Op in Eq. (JT|) being from the VNC 
or from the ANC. For E xx (q,r) these individual contri¬ 
butions, along with their sum, are displayed separately. 
Both E xx {q 1 t) and E xy (q, r) response functions obtained 
with one-body terms only in the NC are substantially in¬ 
creased when two-body terms are also retained. This 
enhancement is found not only at low r, thus corrobo¬ 
rating the sum-rule predictions of Ref. |J], but in fact 
extends over the whole r region studied here. Moreover, 
in the case of the transverse response it affects, in rela¬ 
tive terms, the individual (VNC-VNC) and (ANC-ANC) 
contributions about equally. 


The VNC consists of a linear combination of the isoscalar 
and isovector components of the electromagnetic cur¬ 
rent, weighted respectively by the factors —2 sin 2 9w 
and (1 — 2 sin 2 f?ty) with 6w being the Weinberg an¬ 
gle. The excess transverse strength induced by two-body 
terms in the VNC is consistent with that found in the 
transverse electromagnetic response, and is confirmed by 
experiment as Fig. |T] demonstrates. The two-body en¬ 
hancement in the (ANC-ANC) contribution of E xx (q 1 r) 
is substantial at these relatively large q' s. It decreases 
significantly (for r > 0.01 MeV -1 ) as q is reduced [21] . 
consistently with what is found in calculations of low 
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FIG. 2. (Color online) Euclidean neutral-weak transverse 
(top panel) and interference (lower panel) response functions 
(a/3 = xx and xy in the notation of Ref. [I]) of 12 C at q = 570 
MeV. See text for further explanations. 


q charge-changing weak transitions to specific low-lying 
states, such as the /3-decays and electron and muon cap¬ 
tures studied in Refs. [25, J2BJ, where it amounts to a 
few percent. In principle, the enhancement in the quasi¬ 
elastic region could be measured in parity-violating in¬ 
clusive (e, e') scattering at backward angles. However, 
the smallness of the factor (1 — 4 sin 2 6w), to which the 
relevant (VEM-ANC) interference response function is 
proportional, makes experiments of this type extremely 
difficult. 

In order to obtain more detailed information on the 
energy dependence of the R a p(q,io) response, we em¬ 
ploy the maximum entropy (MaxEnt) method to invert 
E a p (q, t). We describe the method here very briefly, sev¬ 
eral standard references are available [151 IB]. The nu¬ 
merical inversion of a Laplace transform E a p(q,r) with 
its associated statistical errors is a notoriously ill-posed 
problem. The fact that we are interested in the (smooth) 
response around the quasi-elastic peak rather than iso¬ 
lated peaks makes it somewhat more practical. The 
MaxEnt method is based on Bayesian statistical infer¬ 
ence: the “most probable” response function is the one 
that maximizes the posterior probability Pr[i?|i3], i.e., 
the conditional probability of R given E. Bayes theo¬ 
rem states that the posterior probability is proportional 


to the product Pr[i3|i?] x Pr[i?], where Pr[E|i?] is the 
likelihood function and Pr[i? ] is the prior probability. Ar¬ 
guments based on the central limit theorem show that 
the asymptotic limit of the likelihood function is given 
by Pr[E|i?] a exp(—y 2 /2) with y 2 defined as follows. 
Let N t and N u be the numbers of grid points in the 
variables r and ui, respectively. Then the Laplace trans¬ 
form in Eq. © reads (the (^-dependence and subscripts 
a/3 of E a fj{q 1 t) and R a p(q, r) are suppressed for simplic¬ 
ity hereafter) 


N w 

Ei = Y^ Kij Rj , (4) 

3=1 

where = exp (—TiUj) and Rj = A ujjR(uij), and the 
X' 2 follows from 

N r __ _ 

X 2 = E - E 0 (<?”% ~ E i) > (5) 

*.i=i 


where the Ei are obtained from Eq. ©, the Ei are the 
GFMC calculated values, and C is the covariance matrix. 
Therefore, maximizing the likelihood function reduces to 
finding a set of Ri values that minimizes the \ 2 . The 
GFMC errors on Ei are strongly correlated in r, as in¬ 
dividual steps involve only small spatial distances and 
evolutions of the spin-isospin amplitudes. It is therefore 
of paramount importance to estimate the covariance ma¬ 
trix C. 

Limiting ourselves only to the y 2 minimization would 
implicitly be making the assumption that the prior prob¬ 
ability is either unimportant or unknown. However, since 
the response function is positive definite and normal¬ 
izable, it can be interpreted as yet another probability 
function. The principle of maximum entropy states that 
the values of a probability function are to be assigned by 
maximizing the entropy 


S = 


JV„ 

E 

i=l 


[-R(wi) - M(cvi) - R(ui)ln[R(u}i)/M(ui)\ 


A u>i 


( 6 ) 

where the positive definite function M(u>) is the default 
model. It is worthwhile mentioning that the above ex¬ 
pression is applicable even when R(co) and M (w) have dif¬ 
ferent normalizations. The entropy measures how much 
the response function differs from the model. It vanishes 
when R{uj) = M(w), and is negative when R(oj) ^ M(uf). 
The maximum entropy method adds to the simple y 2 
minimization the use of the prior information that the 
response function can be interpreted as a probability dis¬ 
tribution function. We employ historic maximum en¬ 
tropy by minimizing aS — x 2 /2 with the parameter a 
adjusted to make the x 2 equal to one. While more re¬ 
fined methods relying on Bayes statistical inference have 
been developed, we found historic maximum entropy to 
be simple to implement and adequate for our purposes. 
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As a first case we consider the electromagnetic response 
of 4 He. We generated a set of Ne — 2500 GFMC esti¬ 
mates of the Euclidean response functions, obtained from 
independent imaginary-time propagations on a grid of r 
points uniformly distributed between 0 to 0.05 MeV -1 
with At = 0.0005 MeV" 1 . Let E\ n) = E^ n \ Ti ) be 
the Euclidean response function corresponding to the n th 
GFMC propagation. The average Euclidean response 
function and covariance matrix elements are given by 


Ne 


C« = 


n= 1 

i 


Ne 


Ne(N e - 1)^ 


Ei - E, 


(n) 


Ej E)“ 


(7) 


( 8 ) 


In general, the covariance matrix is non-diagonal because 
of correlations between different r,, and the full expres¬ 
sion for the x 2 in Eq. ([ 5 ]) has been used. 

The 4 He longitudinal and transverse response func¬ 
tions (at q = 600 MeV), obtained from inversion of 
El(q,t) and -Et(<Z, t), are shown in Fig. [3j The inver¬ 
sions are, to a very large degree, insensitive to the choice 
of default model response [24] . Results obtained with 
one-body only (dashed line) and (one+two)-body (solid 
line) currents are compared with experimental world 
data [TO] (empty circles). There is excellent agreement 
between the full theory and experiment. Two-body cur¬ 
rents significantly enhance the transverse response func¬ 
tion, not only in the dip region, but also in the quasi¬ 
elastic peak and threshold regions, providing the missing 
strength needed to reproduce the experimental results. 
The band in Fig. [3] provides an estimate for the depen¬ 
dence of the full results on the adopted default model, 
either a flat M{u>) or a gaussian one [24]. The model 
dependence is quite small. 

On the basis of the present 4 He and 12 C calculations, 
a consistent picture of the electroweak response of nu¬ 
clei emerges, in which two-body terms in the nuclear 
electroweak current are seen to produce significant ex¬ 
cess transverse strength from threshold to the dip re¬ 
gion and beyond. Such a picture is at variance with 
the conventional one of inclusive quasi-elastic scatter¬ 
ing, in which single-nucleon knockout is expected to 
be the dominant process in this regime. With the ex¬ 
ception of the leading relativistic corrections contained 
in the nuclear electroweak currents (see Ref. |T]), the 
present calculations are based on a nonrelativistic ap¬ 
proach. Naive kinematical considerations would lead one 
to expect the quasi-elastic peak position in Fig. [3] to be 
at q 2 / (2 m) + A E ~ 211 MeV for q = 600 MeV—we take 
A E ~ 20 MeV to be the separation energy of 4 He into 
a 3+1 cluster. The calculated response functions appear 
to peak at lower w, in fact close to w qe +A E ~ 195 MeV. 
The width of the quasi-elastic peak is also seen to be 


correctly reproduced—the nonrelativistic Fermi gas fails 
to predict this quantity at momentum transfers q ~ 600 
MeV as in Fig. [3] Thus, even at these relatively high 
momentum and energy transfers, the nonrelativistic dy¬ 
namical framework adopted here may be more robust 
than comparisons between nonrelativistic and relativis¬ 
tic Fermi gas models would lead one to conclude EH- 



FIG. 3. (Color online) Electromagnetic longitudinal (top 
panel) and transverse (lower panel) response functions of 4 He 
at q = 600 MeV. Experimental data are from Ref. m- 


A direct evaluation of the 12 C response functions via 
these same methods would require about 100 million core 
hours. We are examining improved methods including 
the use of correlated sampling that could improve the ef¬ 
ficiency of this inversion. We are also exploring methods 
to extend these results to larger nuclei. 
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